/*  -- translated by f2c (version 20100827).
   You must link the resulting object file with libf2c:
	on Microsoft Windows system, link with libf2c.lib;
	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
	or, if you install libf2c.a in a standard place, with -lf2c -lm
	-- in that order, at the end of the command line, as in
		cc *.o -lf2c -lm
	Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,

		http://www.netlib.org/f2c/libf2c.zip
*/

#include "f2c.h"

/* Subroutine */ int igraphdlarrb_(integer *n, doublereal *d__, doublereal *lld, 
	integer *ifirst, integer *ilast, doublereal *rtol1, doublereal *rtol2,
	 integer *offset, doublereal *w, doublereal *wgap, doublereal *werr, 
	doublereal *work, integer *iwork, doublereal *pivmin, doublereal *
	spdiam, integer *twist, integer *info)
{
    /* System generated locals */
    integer i__1;
    doublereal d__1, d__2;

    /* Builtin functions */
    double log(doublereal);

    /* Local variables */
    integer i__, k, r__, i1, ii, ip;
    doublereal gap, mid, tmp, back, lgap, rgap, left;
    integer iter, nint, prev, next;
    doublereal cvrgd, right, width;
    extern integer igraphdlaneg_(integer *, doublereal *, doublereal *, doublereal *
	    , doublereal *, integer *);
    integer negcnt;
    doublereal mnwdth;
    integer olnint, maxitr;


/*  -- LAPACK auxiliary routine (version 3.2) --   
    -- LAPACK is a software package provided by Univ. of Tennessee,    --   
    -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--   
       November 2006   


    Purpose   
    =======   

    Given the relatively robust representation(RRR) L D L^T, DLARRB   
    does "limited" bisection to refine the eigenvalues of L D L^T,   
    W( IFIRST-OFFSET ) through W( ILAST-OFFSET ), to more accuracy. Initial   
    guesses for these eigenvalues are input in W, the corresponding estimate   
    of the error in these guesses and their gaps are input in WERR   
    and WGAP, respectively. During bisection, intervals   
    [left, right] are maintained by storing their mid-points and   
    semi-widths in the arrays W and WERR respectively.   

    Arguments   
    =========   

    N       (input) INTEGER   
            The order of the matrix.   

    D       (input) DOUBLE PRECISION array, dimension (N)   
            The N diagonal elements of the diagonal matrix D.   

    LLD     (input) DOUBLE PRECISION array, dimension (N-1)   
            The (N-1) elements L(i)*L(i)*D(i).   

    IFIRST  (input) INTEGER   
            The index of the first eigenvalue to be computed.   

    ILAST   (input) INTEGER   
            The index of the last eigenvalue to be computed.   

    RTOL1   (input) DOUBLE PRECISION   
    RTOL2   (input) DOUBLE PRECISION   
            Tolerance for the convergence of the bisection intervals.   
            An interval [LEFT,RIGHT] has converged if   
            RIGHT-LEFT.LT.MAX( RTOL1*GAP, RTOL2*MAX(|LEFT|,|RIGHT|) )   
            where GAP is the (estimated) distance to the nearest   
            eigenvalue.   

    OFFSET  (input) INTEGER   
            Offset for the arrays W, WGAP and WERR, i.e., the IFIRST-OFFSET   
            through ILAST-OFFSET elements of these arrays are to be used.   

    W       (input/output) DOUBLE PRECISION array, dimension (N)   
            On input, W( IFIRST-OFFSET ) through W( ILAST-OFFSET ) are   
            estimates of the eigenvalues of L D L^T indexed IFIRST throug   
            ILAST.   
            On output, these estimates are refined.   

    WGAP    (input/output) DOUBLE PRECISION array, dimension (N-1)   
            On input, the (estimated) gaps between consecutive   
            eigenvalues of L D L^T, i.e., WGAP(I-OFFSET) is the gap between   
            eigenvalues I and I+1. Note that if IFIRST.EQ.ILAST   
            then WGAP(IFIRST-OFFSET) must be set to ZERO.   
            On output, these gaps are refined.   

    WERR    (input/output) DOUBLE PRECISION array, dimension (N)   
            On input, WERR( IFIRST-OFFSET ) through WERR( ILAST-OFFSET ) are   
            the errors in the estimates of the corresponding elements in W.   
            On output, these errors are refined.   

    WORK    (workspace) DOUBLE PRECISION array, dimension (2*N)   
            Workspace.   

    IWORK   (workspace) INTEGER array, dimension (2*N)   
            Workspace.   

    PIVMIN  (input) DOUBLE PRECISION   
            The minimum pivot in the Sturm sequence.   

    SPDIAM  (input) DOUBLE PRECISION   
            The spectral diameter of the matrix.   

    TWIST   (input) INTEGER   
            The twist index for the twisted factorization that is used   
            for the negcount.   
            TWIST = N: Compute negcount from L D L^T - LAMBDA I = L+ D+ L+^T   
            TWIST = 1: Compute negcount from L D L^T - LAMBDA I = U- D- U-^T   
            TWIST = R: Compute negcount from L D L^T - LAMBDA I = N(r) D(r) N(r)   

    INFO    (output) INTEGER   
            Error flag.   

    Further Details   
    ===============   

    Based on contributions by   
       Beresford Parlett, University of California, Berkeley, USA   
       Jim Demmel, University of California, Berkeley, USA   
       Inderjit Dhillon, University of Texas, Austin, USA   
       Osni Marques, LBNL/NERSC, USA   
       Christof Voemel, University of California, Berkeley, USA   

    =====================================================================   



       Parameter adjustments */
    --iwork;
    --work;
    --werr;
    --wgap;
    --w;
    --lld;
    --d__;

    /* Function Body */
    *info = 0;

    maxitr = (integer) ((log(*spdiam + *pivmin) - log(*pivmin)) / log(2.)) + 
	    2;
    mnwdth = *pivmin * 2.;

    r__ = *twist;
    if (r__ < 1 || r__ > *n) {
	r__ = *n;
    }

/*     Initialize unconverged intervals in [ WORK(2*I-1), WORK(2*I) ].   
       The Sturm Count, Count( WORK(2*I-1) ) is arranged to be I-1, while   
       Count( WORK(2*I) ) is stored in IWORK( 2*I ). The integer IWORK( 2*I-1 )   
       for an unconverged interval is set to the index of the next unconverged   
       interval, and is -1 or 0 for a converged interval. Thus a linked   
       list of unconverged intervals is set up. */

    i1 = *ifirst;
/*     The number of unconverged intervals */
    nint = 0;
/*     The last unconverged interval found */
    prev = 0;
    rgap = wgap[i1 - *offset];
    i__1 = *ilast;
    for (i__ = i1; i__ <= i__1; ++i__) {
	k = i__ << 1;
	ii = i__ - *offset;
	left = w[ii] - werr[ii];
	right = w[ii] + werr[ii];
	lgap = rgap;
	rgap = wgap[ii];
	gap = min(lgap,rgap);
/*        Make sure that [LEFT,RIGHT] contains the desired eigenvalue   
          Compute negcount from dstqds facto L+D+L+^T = L D L^T - LEFT   

          Do while( NEGCNT(LEFT).GT.I-1 ) */

	back = werr[ii];
L20:
	negcnt = igraphdlaneg_(n, &d__[1], &lld[1], &left, pivmin, &r__);
	if (negcnt > i__ - 1) {
	    left -= back;
	    back *= 2.;
	    goto L20;
	}

/*        Do while( NEGCNT(RIGHT).LT.I )   
          Compute negcount from dstqds facto L+D+L+^T = L D L^T - RIGHT */

	back = werr[ii];
L50:
	negcnt = igraphdlaneg_(n, &d__[1], &lld[1], &right, pivmin, &r__);
	if (negcnt < i__) {
	    right += back;
	    back *= 2.;
	    goto L50;
	}
	width = (d__1 = left - right, abs(d__1)) * .5;
/* Computing MAX */
	d__1 = abs(left), d__2 = abs(right);
	tmp = max(d__1,d__2);
/* Computing MAX */
	d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp;
	cvrgd = max(d__1,d__2);
	if (width <= cvrgd || width <= mnwdth) {
/*           This interval has already converged and does not need refinement.   
             (Note that the gaps might change through refining the   
              eigenvalues, however, they can only get bigger.)   
             Remove it from the list. */
	    iwork[k - 1] = -1;
/*           Make sure that I1 always points to the first unconverged interval */
	    if (i__ == i1 && i__ < *ilast) {
		i1 = i__ + 1;
	    }
	    if (prev >= i1 && i__ <= *ilast) {
		iwork[(prev << 1) - 1] = i__ + 1;
	    }
	} else {
/*           unconverged interval found */
	    prev = i__;
	    ++nint;
	    iwork[k - 1] = i__ + 1;
	    iwork[k] = negcnt;
	}
	work[k - 1] = left;
	work[k] = right;
/* L75: */
    }

/*     Do while( NINT.GT.0 ), i.e. there are still unconverged intervals   
       and while (ITER.LT.MAXITR) */

    iter = 0;
L80:
    prev = i1 - 1;
    i__ = i1;
    olnint = nint;
    i__1 = olnint;
    for (ip = 1; ip <= i__1; ++ip) {
	k = i__ << 1;
	ii = i__ - *offset;
	rgap = wgap[ii];
	lgap = rgap;
	if (ii > 1) {
	    lgap = wgap[ii - 1];
	}
	gap = min(lgap,rgap);
	next = iwork[k - 1];
	left = work[k - 1];
	right = work[k];
	mid = (left + right) * .5;
/*        semiwidth of interval */
	width = right - mid;
/* Computing MAX */
	d__1 = abs(left), d__2 = abs(right);
	tmp = max(d__1,d__2);
/* Computing MAX */
	d__1 = *rtol1 * gap, d__2 = *rtol2 * tmp;
	cvrgd = max(d__1,d__2);
	if (width <= cvrgd || width <= mnwdth || iter == maxitr) {
/*           reduce number of unconverged intervals */
	    --nint;
/*           Mark interval as converged. */
	    iwork[k - 1] = 0;
	    if (i1 == i__) {
		i1 = next;
	    } else {
/*              Prev holds the last unconverged interval previously examined */
		if (prev >= i1) {
		    iwork[(prev << 1) - 1] = next;
		}
	    }
	    i__ = next;
	    goto L100;
	}
	prev = i__;

/*        Perform one bisection step */

	negcnt = igraphdlaneg_(n, &d__[1], &lld[1], &mid, pivmin, &r__);
	if (negcnt <= i__ - 1) {
	    work[k - 1] = mid;
	} else {
	    work[k] = mid;
	}
	i__ = next;
L100:
	;
    }
    ++iter;
/*     do another loop if there are still unconverged intervals   
       However, in the last iteration, all intervals are accepted   
       since this is the best we can do. */
    if (nint > 0 && iter <= maxitr) {
	goto L80;
    }


/*     At this point, all the intervals have converged */
    i__1 = *ilast;
    for (i__ = *ifirst; i__ <= i__1; ++i__) {
	k = i__ << 1;
	ii = i__ - *offset;
/*        All intervals marked by '0' have been refined. */
	if (iwork[k - 1] == 0) {
	    w[ii] = (work[k - 1] + work[k]) * .5;
	    werr[ii] = work[k] - w[ii];
	}
/* L110: */
    }

    i__1 = *ilast;
    for (i__ = *ifirst + 1; i__ <= i__1; ++i__) {
	k = i__ << 1;
	ii = i__ - *offset;
/* Computing MAX */
	d__1 = 0., d__2 = w[ii] - werr[ii] - w[ii - 1] - werr[ii - 1];
	wgap[ii - 1] = max(d__1,d__2);
/* L111: */
    }
    return 0;

/*     End of DLARRB */

} /* igraphdlarrb_ */

